home *** CD-ROM | disk | FTP | other *** search
/ AOL File Library: 2,801 to 2,900 / aol-file-protocol-4400-2801-to-2900.zip / AOLDLs / C++ Files Library / FFT and Complex Numbers Library / FFTLibrary.sit / FFT library ƒ / IFFT.cp < prev    next >
Text File  |  1995-01-15  |  3KB  |  96 lines
  1. //
  2. //    File: IFFT.cp
  3. //    Copyright ⌐1994-1995 James Wilson
  4. //    All rights reserved.
  5. //
  6. //    Revisions:
  7. //        01-15-95: Removed the line that zeros the imaginary during the
  8. //        inverse transformation. This only slows the code, and does not
  9. //        take in the possibility of complex data. If the data is real,
  10. //        code that uses this library should ignore the imaginary results
  11. //        of the IFFT.
  12. //        08-02-94: Took out the constant multiplication part of arg and
  13. //        placed into halfArg. halfArg is multiplied by p to get arg.
  14. //        Removed the "f" suffix from 1.0 on line 34. Removed the "f"
  15. //        suffix from 0.0 on line 86. Replaced the name "temp" by "swap."
  16. //        07-29-94: Replaced all instances of floating-point division with
  17. //        multiplication by reciprocals, since multiplication is much
  18. //        faster on PowerPC.
  19. //        
  20. //    This is the implementation file for the Inverse Fourier transform (IFFT)
  21. //    algorithm. The algorithm is the same as the FFT, except for the sign
  22. //    of the comlex exponential and the final multiplication by 1/n. I assume
  23. //    that the reader or user of this code has a decent understanding of both
  24. //    the DFT, FFT, and inverse transforms.
  25.  
  26. #ifndef __FFT__
  27. #include "FFT.h"
  28. #endif
  29.  
  30. void IFFT(const long n, const long nu, ComplexArray functionOut)
  31. {
  32.     long                    n2, nu1, i, k, p, x, tempIndex, denom;
  33.     double                    arg, reciprocalN, halfArg;
  34.     Complex                    exponent, swap;
  35.     
  36.     // Initalization process.
  37.     x = 1L;
  38.     n2 = n / 2L;
  39.     nu1 = nu - 1L;
  40.     k = 0L;
  41.     reciprocalN = 1.0 / n;
  42.     halfArg = twoPi * reciprocalN;
  43.     
  44.     // First step: calculate the Fourier transforms of the ComplexArray values.
  45.     // The results are scrambled and are unscrambled at the end (see
  46.     // below).
  47.     while(x <= nu)
  48.     {
  49.         for(i = 1L; i <= n2; i++, k++)
  50.         {
  51.             denom = LongPower(2L, nu1);
  52.             p = BitReverse(k/denom, nu);
  53.             arg = halfArg * p;
  54.             
  55.             exponent.Real() = cos(arg);
  56.             exponent.Imag() = sin(arg);
  57.             
  58.             tempIndex = k + n2;
  59.             
  60.             swap = exponent * functionOut[tempIndex];
  61.             functionOut[tempIndex] = functionOut[k] - swap;
  62.             functionOut[k] += swap;
  63.         }
  64.         
  65.         k += n2;
  66.         
  67.         if(k >= n)
  68.         {
  69.             x++;
  70.             n2 /= 2L;
  71.             nu1--;
  72.             k = 0L;
  73.         }
  74.     }
  75.     
  76.     // Last step: unscramble the results. This is done by swapping the
  77.     // values in functionOut[k] with functionOut[BitReverse(k, nu)].
  78.     // Then multiply each element by 1/n. The imaginary part will be close
  79.     // to zero, but due to precision problems, it will not be equal to 0.0.
  80.     while(k < n)
  81.     {
  82.         i = BitReverse(k, nu);
  83.         
  84.         if(i > k)
  85.         {
  86.             swap = functionOut[k];
  87.             functionOut[k] = functionOut[i];
  88.             functionOut[i] = swap;
  89.         }
  90.         
  91.         functionOut[k].Real() *= reciprocalN;
  92.         
  93.         k++;
  94.     }
  95. }
  96.